/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#include <math.h>
#include <string.h>
#include <stdint.h>

#include "fastonebigheader.h"

#define log_polytomous_rasch(x, beta, taufunc, n_tau, dbeta, dtau)      \
  ({                                                                    \
     float logfactors[1 + n_tau];                                       \
     float dlogfactors_dbeta[1 + n_tau];                                \
     float dlogfactors_dtau[1 + n_tau][n_tau];                          \
     unsigned int argmax;                                               \
                                                                        \
     logfactors[0] = 0;                                                 \
     dlogfactors_dbeta[0] = 0;                                          \
     memset (dlogfactors_dtau[0], 0, (n_tau) * sizeof (float));         \
     argmax = 0;                                                        \
                                                                        \
     for (unsigned int _lpri = 0; _lpri < (n_tau); ++_lpri)             \
       {                                                                \
         logfactors[_lpri + 1] =                                        \
           logfactors[_lpri] + (beta - taufunc (_lpri));                \
         dlogfactors_dbeta[_lpri + 1] =                                 \
           dlogfactors_dbeta[_lpri] + 1;                                \
         memcpy (dlogfactors_dtau[_lpri + 1],                           \
                 dlogfactors_dtau[_lpri],                               \
                 (n_tau) * sizeof (float));                             \
         dlogfactors_dtau[_lpri + 1][_lpri] -= 1;                       \
                                                                        \
         if (logfactors[_lpri + 1] > logfactors[argmax])                \
           {                                                            \
             argmax = _lpri + 1;                                        \
           }                                                            \
       }                                                                \
                                                                        \
     /*
      *    exp (logfactors[x])          exp (logfactors[x] - logfactors[argmax])
      * ------------------------- = -----------------------------------------------
      * sum_y exp (logfactors[y])   sum_y exp (logfactors[y] - logfactors[argmax])
      */                                                                \
                                                                        \
     float denominator = 0;                                             \
     float ddenominator_dbeta = 0;                                      \
     float ddenominator_dtau[n_tau];                                    \
                                                                        \
     memset (ddenominator_dtau, 0, (n_tau) * sizeof (float));           \
                                                                        \
     for (unsigned int _lpri = 0; _lpri < 1 + (n_tau); ++_lpri)         \
       {                                                                \
         float _lprdtd =                                                \
           fastexp (logfactors[_lpri] - logfactors[argmax]);            \
         denominator += _lprdtd;                                        \
         ddenominator_dbeta += _lprdtd *                                \
           (dlogfactors_dbeta[_lpri] - dlogfactors_dbeta[argmax]);      \
                                                                        \
         for (unsigned int _lprdtk = 0; _lprdtk < (n_tau); ++_lprdtk)   \
           {                                                            \
             ddenominator_dtau[_lprdtk] += _lprdtd *                    \
               (dlogfactors_dtau[_lpri][_lprdtk] -                      \
                dlogfactors_dtau[argmax][_lprdtk]);                     \
           }                                                            \
       }                                                                \
                                                                        \
     for (unsigned int _lprdtk = 0; _lprdtk < (n_tau); ++_lprdtk)       \
       {                                                                \
         dtau[_lprdtk] = dlogfactors_dtau[x][_lprdtk] -                 \
                         dlogfactors_dtau[argmax][_lprdtk] -            \
                         ddenominator_dtau[_lprdtk] / denominator;      \
       }                                                                \
                                                                        \
     dbeta = dlogfactors_dbeta[x] - dlogfactors_dbeta[argmax] -         \
       ddenominator_dbeta / denominator;                                \
                                                                        \
     logfactors[x] - logfactors[argmax] - fastlog (denominator);        \
   })

#define dlog_polytomous_rasch_dall(dtau, x, beta, taufunc, n_tau)       \
  ({                                                                    \
     float logfactors[1 + n_tau];                                       \
     float dlogfactors_dbeta[1 + n_tau];                                \
     float dlogfactors_dtau[1 + n_tau][n_tau];                          \
     unsigned int argmax;                                               \
                                                                        \
     logfactors[0] = 0;                                                 \
     dlogfactors_dbeta[0] = 0;                                          \
     memset (dlogfactors_dtau[0], 0, (n_tau) * sizeof (float));         \
     argmax = 0;                                                        \
                                                                        \
     for (unsigned int _dlprdti = 0; _dlprdti < n_tau; ++_dlprdti)      \
       {                                                                \
         logfactors[_dlprdti + 1] =                                     \
           logfactors[_dlprdti] + (beta - taufunc (_dlprdti));          \
         dlogfactors_dbeta[_dlprdti + 1] =                              \
           dlogfactors_dbeta[_dlprdti] + 1;                             \
         memcpy (dlogfactors_dtau[_dlprdti + 1],                        \
                 dlogfactors_dtau[_dlprdti],                            \
                 (n_tau) * sizeof (float));                             \
         dlogfactors_dtau[_dlprdti + 1][_dlprdti] -= 1;                 \
                                                                        \
         if (logfactors[_dlprdti + 1] > logfactors[argmax])             \
           {                                                            \
             argmax = _dlprdti + 1;                                     \
           }                                                            \
       }                                                                \
                                                                        \
     float denominator = 0;                                             \
     float ddenominator_dbeta = 0;                                      \
     float ddenominator_dtau[n_tau];                                    \
                                                                        \
     memset (ddenominator_dtau, 0, (n_tau) * sizeof (float));           \
                                                                        \
     for (unsigned int _dlprdti = 0; _dlprdti < 1 + n_tau; ++_dlprdti)  \
       {                                                                \
         float _dlprdtd =                                               \
           fastexp (logfactors[_dlprdti] - logfactors[argmax]);         \
         denominator += _dlprdtd;                                       \
         ddenominator_dbeta += _dlprdtd *                               \
           (dlogfactors_dbeta[_dlprdti] - dlogfactors_dbeta[argmax]);   \
                                                                        \
         for (unsigned int _dlrpdtk = 0; _dlrpdtk < n_tau; ++_dlrpdtk)  \
           {                                                            \
             ddenominator_dtau[_dlrpdtk] += _dlprdtd *                  \
               (dlogfactors_dtau[_dlprdti][_dlrpdtk] -                  \
                dlogfactors_dtau[argmax][_dlrpdtk]);                    \
           }                                                            \
       }                                                                \
                                                                        \
     for (unsigned int _dlrpdtk = 0; _dlrpdtk < n_tau; ++_dlrpdtk)      \
       {                                                                \
         dtau[_dlrpdtk] = dlogfactors_dtau[x][_dlrpdtk] -               \
                          dlogfactors_dtau[argmax][_dlrpdtk] -          \
                          ddenominator_dtau[_dlrpdtk] / denominator;    \
       }                                                                \
                                                                        \
     dlogfactors_dbeta[x] - dlogfactors_dbeta[argmax] -                 \
       ddenominator_dbeta / denominator;                                \
   })
